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Abstract 

In this work the path integral formulation for rigid rotors, proposed by Miiser and Berne [Phys. 
Rev. Lett. 77, 2638 (1996)], is described in detail. It is shown how this formulation can be used to 
perform Monte Carlo simulations of water. Our numerical results show that whereas some proper- 
ties of water can be accurately reproduced using classical simulations with an empirical potential 
which, implicitly, includes quantum effects, other properties can only be described quantitatively 
when quantum effects are explicitly incorporated. In particular, quantum effects are extremely 
relevant when it comes to describing the equation of state of the ice phases at low temperatures, 
the structure of the ices at low temperatures, and the heat capacity of both liquid water and the 
ice phases. They also play a minor role in the relative stability of the ice phases. 
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I. INTRODUCTION 



In 1948 Feynman proposed the path integral formulation of quantum mechanicsi'^. In the 
late seventies the works of Barker and of Chandler and Wolynes showed that this formulation 
could be implemented in the statistical mechanical study of condensed matter by performing 
classical simulations of a modified Hamiltonian^'^. It was demonstrated that the partition 
function of a quantum system of particles in its discretised form is formally identical to 
that of a classical system consisting of ring polymers. Thus a number of the techniques 
and methods that had already been derived for classical simulations could now be adapted 
to perform quantum simulations. Since then path integral simulations have been used to 
study the behaviour of a large number of systems. A detailed description of the path integral 
technique in statistical mechanics, and its applications, can be found in several reviews^*^. 

The path integral formulation can be implemented both in Monte Carlo (MC) and in 
molecular dynamics (MD) simulations. MC simulations only provide thermodynamic prop- 
erties. Although there has been some controversy concerning the description of the correct 
dynamics, there are now MD methods, such as centroid molecular dynamics^'^ or the ring 
polymer molecular dynamics^ that approximate the quantum dynamics of the system. 

The path integral method has also been extended to study rigid rotors. The first quantum 
MC simulations of rigid bodies were performed in the mid-eighties using a semi-classical 
approximation to derive the rotational contribution to the partition function^^iii. This 
semi-classical approximation was also used within molecular dynamics in the beginning of 
the ninetiesi^. A few years later Miiser and Berne derived a propagator for the rotational 
contribution to the density function for rigid bodies^*^. Centroid molecular dynamics has 
also been recently extended to deal with rigid rotors^^. These attempts to extend the 
path integral method to rigid bodies are motivated by the desire to describe the quantum 
behaviour of molecules in a computationally efficient way, since various condensed matter 
properties are more likely to be affected by inter-molecular vibrations rather than intra- 
molecular contributions- . In this situation, it seems reasonable to describe the molecules 
as rigid rotors and ignore the intra-molecular vibrations. Usually intra-molecular vibrations 
exhibit high frequencies that require to use a large number of 'replicas' of the system, thus 
increasing considerably the computational cost of quantum simulations. This problem can 
be alleviated by the use of smart techniques such as the recently proposed ring polymer 
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contractor methodic. That said, it is our opinion that the description of molecules as rigid 
bodies is particularly interesting in its self; since it allows one to separate inter-molecular 
and intra-molecular quantum effects, thus assigning their relative influence on the different 
properties of the condensed phases. 

In this work we shall describe the path integral formulation for rigid rotors proposed by 



Miiser and Berne. Even though this formulation has already described in Refs. [13| and [14 , 
the implementation of a path integral Monte Carlo involves many technical details which, 
in our opinion, are worth describing in detail. In addition, we show how this method can be 
applied to study quantum effects for a range of properties of water. 

II. PATH INTEGRAL FOR RIGID ROTORS 
A. General path integral formulation 

The behaviour of a system of quantum non-spherical rigid particles can be described 
by the Schrodinger equation: 

= ^A^A (1) 

where the Hamiltonian is given by: 

H = f *™ + f + U (2) 

where T*'''^ is the translational kinetic energy operator of the centre of mass of the molecules, 
T^°* is the orientational kinetic energy operator of the rigid molecules, and U is the potential 
energy operator. The partition function for this system at inverse temperature /3 can be 
written as: 

Z = 5^(*,|e-^^|^A) = Tr{e-^^) = ^e"^^^ (3) 

A A 

where \E'a and Ex are the eigenf unctions and the eigenvalues, respectively, of the Hamilto- 
nian. The solutions of the Schrodinger equation form a complete orthonormal basis (i.e. 
(\I'a|^a') = f^AA')- The operator p = e~^^ is defined by its Taylor expansion: 

°0 TTk 

k=0 

The second term in Eq. [3]implies first an integration over the coordinates of the system (r, u) 
followed by a sum over the states A. Here r represents the set of Cartesian coordinates of the 



centre of mass of each particle and u are the set of Euler angles that define the orientation 
of each particle in the system. Note that although we have used the eigenfunctions of the 
Hamiltonian to write the partition function, any complete basis set can be used due to the 
fact that the trace is invariant with respect to a change of basis^. 

We now introduce the density function, which in the coordinate space is defined as: 

p(r, uj, r', uj\ P) = J2 ^A(r, co) expi-(3Ex)^ xir' , u') (5) 

A 

The second term in Eq. [3] can be re-written as: 

Z = Tr{e-^^) = l^'$l{r,u)exp{-/3Ex)'^x{r,u)drduj= f p{r,u,r,uj, P)drduj (6) 
J X 

where one first sums over the states A to obtain p(r, u, r, u, /3) and then integrates over the 
coordinates. One interesting property of the density function is that the product of two 
density functions results in another density function: 

j p{T, cu, r", oj"; /3i)p(r", cu", r', co'; (32)dv"du" = p(r, r', uj'- (3, + /^a) (7) 

Naturally this can be generalised to any given number of terms. Using this convolution 
property of the density function the partition function can be re-written as: 

Z = j P {v\ uj\ v\ u?- ^ ...p cu^-i, r V; ^ p (^r^, u;^, r^, u^- ^ dvHuj\ . .dx^ du'' 

(8) 

Usually it is not possible to solve the Schrodinger equation for a system of rigid non- 
spherical interacting particles, so the eigenfunctions \1/a and eigenvalues Ex are unknown. 
However, as mentioned before, the partition function (or what is the same, the trace of the 
density matrix operator p) is invariant independently of the complete set of basis functions 
used. A convenient complete basis set is the one formed by the eigenfunctions of position 
and orientation operators which, for non-spherical particles, is given by: 

|rV) = (5(r - - (9) 

Therefore, we can write the partition function using the eigenfunctions of the position and 
orientation operator by substituting the sum over the states A by an integration over r and 
UJ. In this way, the partition function can be calculated by evaluating the density function 
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p*'*+i(/3/P): 



p'''+\f3/P) = p{r\ uj\ r'+\ u'+\ (3/P) = ( r'u' exp -/3(T*™ + T'"' + U)/P 



J , t i+l , ,t+l 



(10) 

where |r*+^a;*^^) are the eigenf unctions of the position operator. In principle, the exponential 
cannot be factorised because the potential and kinetic energy do not commute. However, 
using the Trotter formular^, one can write the following exact formula in the limit of infinite 
P: 



lim p(r*,a;*,r*+\w*+\/3/P) 

P— )-oo 



(11) 



exp 



-I3U/2P 



exp 



„t+l, ,t+l 

r uj 



where P is the Trotter number. Given the ring polymer isomorphism, P is also known as 
as the number of replicas or beads. 

The operator U is diagonal in the coordinate representation, and for rigid rotors T**"" and 
jirot commute. Therefore, for rigid rotors p*'*"^^(/3/P) can be approximated as: 



p*'*+^(/^/p) ^ p'^:[\p/p)p'^::\(5/p)pror{NP) 



t,t+i. 



t,t+i. 



(12) 



where 



Ppit^iNP) = exp 



exp ( -f3f''yP 



pYot\(ilP) = W exp{-/3T°'/P 



.t+i 



(13) 

(14) 
(15) 



where f/(r*,a;*) is the potential energy of rotors whose positions and orientations are 
specified by r*, w*. 

It can be shown that the translational contribution is given by^: 



p'^^^'iP/P) 



MP 

2TTh^l3 



3/2 



exp 



N 



MP 



2 n 



(16) 



where M is the rotor mass and r* are the coordinates of the centre of mass of replica t of 
rotor i. 

B. Rotational propagator for a free rotor 
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In order to evaluate the partition function for a system of N free rotors it is necessary to 
evaluate the rotational propagator: 



UJ 



exp 



' pTrot 



t+1 



which can be written as the product: 



exp 



' P 



UJ 



t+l 



N 

n 

i=l 



P rinrot 



exp ( --T, 



CO. 



t+l 



(17) 



(18) 



where louf) are the eigenfunctions of the orientation operator, which, as mentioned before, 
are Dirac delta functions: 

\ool) = S{n - n*) (19) 

Here fl represents the three Euler angles (^,0, x). To simplify the notation, we shall drop 
the subindex i, and from here on focus our attention on a single free rotor. 

The eigenfunctions of the angular position can be expanded in a basis set of the 
eigenfunctions of the asymmetric top \JMK) (the derivation of the eigenfunctions of the 
asymmetric top are given in Appendix A): 



J M K 

Using this expansion, the rotational propagator can be rewritten as: 



(20) 



UJ 



exp 



t+l 



JMKj 



(21) 



As the functions | JMK) are the eigenfunctions of the Schrodinger equation for the asym- 
metric top, it follows that: 



P nnrot 



JMK ) = exp ( 



Using this, the rotational propagator can be written as: 



id 



J M K 



P K J 



JMKj 



JMK 



> 



Reordering this expression one has: 



UJ 



exp ( -^T-* 



UJ 



t+l 



(22) 



(23) 



J M K 



(24) 
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The location of the laboratory frame defining the Euler angles is arbitrary. For conve- 
nience we choose a laboratory frame such that the Euler angles of replica t are all zero (i.e., 
= (0, 0, 0)). This change leads to: 



(25) 



and 



= 5(n - (26) 

A tilde is added to Q (i.e., Q) in order to remind ourselves that the Euler angles are defined 
in a laboratory frame in which the Euler angles of bead t are zero. Thus 17*+^ are the Euler 
angles of bead t + 1 in this arbitrary frame. 

To simplify this expression further, the eigenfunctions of the asymmetric top | JMK) are 
expanded in a basis set formed by the eigenfunctions of the symmetric top ( | JMK) ) : 



exp ( --T 



'") = EEE 

J M K 



\JMk) = Y,Af^\JMK) 

K 

Eq. [23] can now be re- written as: 

K 

J M K \ K J ^ ^ 



(27) 



(28) 



K 



The eigenfunctions of the symmetric top are (see Ref. 1191 ) 

'2J + n^/^ 



^JA/i^(^^) 



exp(zM0)d](^^(^) exp(zirx) 



(29) 



where d\.ij^{ff) are the Wigner functions (given in Appendix B). Using this expression we 
have: 

^JMi.(0) = ( ) 4/x(0) = ( ) ^MK (30) 



and 



2J+ 1 



1/2 



Stt^ 



exp(-zM0*+i)rf5[,^(^~*+i) exp(-zirx 



(31) 



having made use of the following relation for Wigner functions: (iMii'(O) = ^mk- Miiser and 
Berne concluded that, substituting these expressions, Eq. [28] could be written as^^: 



exp 



J M K 



p 

2J+1 



UJ 



t+1 



A 



(JM) 
KM 



The propagator is a real quantity and this can be seen by symmetrising it with respect 
to M. This can be achieved by calculating the average of the M and —M contributions for 
any given J and K: 



1 
2 



+ 



2J+ 1 
87r2 
2J+ 1 



A 



(JM) 
'KM 



4Er')diMio'^') 



exp ( --pi^^ J ctMMlf^" ' 1 exp(-iM(0*+^ + x + )) 



A 



(JM) 
KM 



exp 



where M denotes —M. Given that (see Ref. l20f): 



(34) 



and 



.(JM) 
^KM 



A 



(JM) 
KM 



exp 



P k 



Eqj33] can be simplified to: 



2J+ 1 



(JM) 
'KM 



exp 



C^A/Af(^t+l) cos(M(0t+i + Xt+l)) 



(35) 



(36) 



(37) 



Using this result, Miiser and Berne obtained the following final expression for the orienta 
tional propagator: 

oo J J 



t 

Pvot 



r— n A A— _ 7 r> T V / 



(38) 



J=0 M=-J 



where /*'*^!. ^ is a function of the relative Euler angles between beads t and t + 1 and is 
given by: 



t,t+i 

i,J,M,k 



2J+1 



A 



(JM) 
^KM 



</M(^r')cos[M(0r^+xr^)] 



t,t+l I 



(39) 



C. Path integral simulations for a system of rigid rotors 
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Once the three contributions to the density function (EqIT2i). the potential (Eq. [13]), the 
translational (Eq. [T6|) and rotational (Eq. [18] and [38]) have been obtained, the partition 
function of a system of rigid rotors can be calculated by substituting them into Eq. [HI 

MP ^ ^ R ^ \ ^ ^ 

2^ E E « - ^n'- J E f n n ^!:'wn 

^P''' i=i t=i t=i ) i=i i=i 

If one uses a pairwise potential then f/* = ^ . X]j>i U (r-, cu*, r* , u;*). One can now see that 
the partition function of a system of quantum particles is isomorphic to that of a system 
of classical ring polymers. Each replica t of molecule i interacts with the replicas with the 
same index t of the remaining particles through the inter-molecular potential U , and with 
replicas t — 1 and t + 1 of the same molecule i through a harmonic potential, whose coupling 
parameter depends on the mass of the molecules (M) and on the temperature (/3 = l/ksT), 
and through terms p*ot~/ pltt'l which incorporate the quantisation of the rotation and 
which depends on the relative orientations of replicas t — 1 and t, and t and t + 1. The 
convention of the Euler angles used in the present calculations, as well as the conversion 
from Cartesian coordinates to Euler angles, are given in Appendix C. The procedure to 
obtain the Euler angles of replica t + 1 in the body frame of replica t is outlined in Appendix 
D. 

Note that the function pHtlW/^) (Eqs. 38 and 39) depends solely on two angles, ^*'*+^ 
and 0*'*+^ + It is convenient to compute the density function for a grid of values over 

the angles 9 and + X ^ind save these data in tabular form prior to any simulation. The 
value of the density function for any particular 9 and + x can then estimated using an 
interpolation algorithm in conjunction with the tabulated data. 

The internal energy E can be calculated from the partition function Z using the thermo- 
dynamic relation: 

By substituting the partition function it can be shown that the internal energy can be 
calculated as: 

E = Ktra + Krot + (42) 



where the functional forms of these three components are: 



rot 




) 



_Ee^JM) 

e p K 




(45) 



(43) 



(44) 



As with the rotational contribution to the density function, the numerator of the last term 
in Eq. |44] was calculated prior to simulation for a grid of angles 9 and (f) + X and saved in 
tabular form. 

The partition function for the NpT ensemble can be calculated using: 



The implementation of the NpT ensemble in PIMC has already been discussed in previous 
works2ii2S. 

III. QUANTUM EFFECTS IN WATER 

In this section we shall see how the path integral formulation for rigid rotors can be used 
to study quantum effects related to the atomic mass in water. Water exhibits such quantum 
effects even at room temperature. For example, properties such as the melting temperature, 
the temperature of maximum density (TMD), or the heat capacity at constant pressure, all 
show changes upon isotopic substitution^"^. Experimental data for different water isotopes 
are given in Table HI In particular, the melting temperature of water is about 4.49K higher 
in tritiated water than in H2O. A more dramatic effect can be seen in the increase in the 
TMD, which is about 8.91K higher in tritiated water. In contrast, isotopic substitution of 
the oxygen mass has little influence on the properties of liquid water (see Table [I]). This 
indicates that quantum effects are mainly due to the light mass of hydrogen, which leads 
to small values of the principal moments of inertia. The lowest principal moment of inertia 
increases when hydrogen mass increases, whereas it remains almost unchanged when oxygen 
is substituted. The effect that this has on the magnitude of quantum effects can be seen 




(46) 
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by using the following approximate expression to estimate the quantum effects of a rigid 
asymmetric top2^: 



N 2A{kBTY 



{F') (rj) (r|) (r^) 

M I A Ib Ic 



24 ^ \Ia IbIc) 



(47) 



cyclic 

where F is the force that acts on the centre of mass of the particles, 1b and Iq are the 
principal moments of inertia of the particles, and F^, V b and Vc are the torques associated 
with each principal axis of inertia. This expression indicates that both the low moments 
of inertia as well as the strength of the hydrogen bond (which leads to high values of the 
average torque on the molecules) are responsible for the importance of quantum effects in 
water. 

The importance of nuclear quantum effects, as well as the limitations of performing 
classical simulations for water, have already been discussed by Stillinger and Rahman in 
1974, one of the first articles published concerning simulations of water^^. In this pioneering 
paper it was pointed out that classical models are likely to overestimate the temperature 
difference between the TMD and the melting temperature of water, which experimentally is 
3.98K. They speculated that this difference would be about 14K for classical models, which 
was obtained by plotting the difference between the TMD and the melting temperature for 
water, deuterated water and tritiated water, as a function of the inverse of the hydrogen mass 
and then linearly interpolating to infinite mass, which would correspond to the classical limit. 
This hypothesis was also reiterated more recently by Guillot in his review of water models^-. 
A thorough investigation of classical rigid non-polarisable models eventually showed that 
the temperature difference between the TMD and the melting temperature was even larger; 
approximately 30K— . 

The first simulations of water to explicitly include quantum effects were performed almost 
thirty years ago^Si"— i^. In these works, quantum effects in liquid water and deuterated water 
were investigated by using the rigid non-polarisable ST2 model. It was observed that liquid 
water is less structured than heavy water, and that classical water (which corresponds to the 
limit of infinite mass) is more structured than either of these isotopes. In the same period, 
Wallqvist and Berne studied quantum effects using a flexible model^. 

Since these seminal works, a good number of publications have appeared treating the 
quantum effects in water, the majority focusing on the liquid phase^^"^"— . Some of these 
studies implemented potentials designed to be used within classical simulations. Even 
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TABLE I. Experimental data for different water isotopes. The heat capacity at constant pressure 
is given at T=290K. 



IH2O 2H2O ^HaO iHa^^O 



T^eit (K) 273.15 276.97 277.64 273.43 

TMD (K) 277.13 284.33 286.55 277.35 

PTAID X 102 (molec./A3) 3.344 3.326 3.322 3.347 

TMD-Tmeit (K) 3.98 7.36 8.91 3.92 

Cp (cal mol-^ K-i) 17.93 19.60 

though this approach is vahd when it comes to studying the effect that inclusion of quantum 
effects have on the diverse properties of water, a quantitative description requires the use of 
a potential that has been specifically designed to use within path integral simulations. So far 
a few water models have been developed expressly to be used within quantum simulations. 
These include both rigid^- and flexible models that have either been fitted to experimen- 
tal data^>^ or to more accurate ab initio calculations^. Recently more sophisticated path 
integral Car- Parr inello molecular simulations of liquid water and ice Ih have also been under- 
taken, which showed that ab initio calculations are also improved by incorporating nuclear 
quantum effects^. As far as we are aware, there are only a few studies that have considered 
quantum effects for ice I/j— i^"— . There are also a couple of studies that have examined the 
liquid-solid^'^ and solid- vapour interface^. 

In this work we treat the water molecule as being rigid. This means that our simulations 
are only capable of providing information about the low frequency inter-molecular librations 
(for water these are below 900cm~^), whereas high frequency intra-molecular vibrations 
(between 1500 cm~^ for bending and 3500 cm~^ for stretching) will be ignored. The number 
of replicas required to accurately reproduce the properties of a quantum system depends on 
the largest vibration frequency of the system: 

P > (48) 

From this expression it can be seen that at room temperature approximately P =30 replicas 
should be used for flexible models (for which the higher frequencies are around 3500 cm~^), 
whereas P =5 or 6 is sufficient for a rigid model (for which the higher frequencies are 
about 1000 cm~^). This permits a considerable reduction in the computational cost of the 
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simulations. 



IV. METHODOLOGY 



In the quantum simulations presented in this work, water was described using the recently 
proposed TIP4PQ/2005 model^^^, which is the quantum counterpart of the TIP4P/2005 
model^. The classical TIP4P/2005 model was found to provide the best overall descrip- 
tion of water from among the many simple rigid non-polarisable models available in the 
literature^"—. In both models a Lennard- Jones (LJ) centre is located on the oxygen site, 
positive charges on the hydrogens and a negative charge along the bisector of the oxygen- 
hydrogen vectors. The total energy of the system is given by: 

^ = EE^ 

I j>i L 

where Vij represents the distance between the oxygen atoms in molecules i and j, Vmn is 
the distance between the charge qm of molecule i and charge g„ of molecule j. a and 
e are the LJ parameters. The parameters of the TIP4P/2005 and TIP4PQ/2005 models 
are given in Table [Tll The only difference between the two models is that the hydrogen 
charges have been increased in the TIP4PQ/2005 model by 0.02e, which leads to stronger 
electrostatic interactions. This compensates for the loss of structure and the increase in 
energy that is observed when quantum simulations are performed. The same recipe was 
also used previously to obtain quantum counterparts of the classical SPC/F— and TIP5P'^- 
models. The increase in the charges enhances the dipole moment of the water molecule 
from 2.305D in TIP4P/2005 to 2.380D in TIP4PQ/2005 (higher multipole moments will 
obviously also change accordingly). Note that in some cases for flexible models the change 
in geometry caused by the incorporation of quantum effects with respect to the classical 
limit also leads to an enhancement of the dipole moment of the water molecule. Therefore, 
for flexible models it is generally not necessary to increase the charges in order to perform 
quantum simulations^^. 

In this work, the influence of quantum effects in water was investigated by perform- 
ing NpT PIMC simulations using the formulation for rigid rotors proposed by Miiser and 
Berne^^ described previously. Classical NpT simulations were also performed for com- 
parison. The simulation box contained 300 water molecules for the liquid phase, and 432 
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EE 



(49) 



TABLE II. Parameters of the models TIP4P/2005 and TIP4PQ/2005. 



Model a{k) 


e/kB (K) Z HOH (deg 


) doH 


doM 


QH (e) 


TIP4P/2005 3.1589 


93.2 


104.52 


0.9572 


0.1546 


0.5564 


TIP4PQ/2005 3.1589 


93.2 


104.52 


0.9572 


0.1546 


0.5764 



molecules for ices Ih and II. The initial proton disordered configuration of ice Ih was obtained 
using the algorithm of Buch et al.—^. The LJ interaction was truncated at 8.5 A. Standard 
long-range corrections for the LJ part of the potential were added. Coulombic interactions 
were calculated using Ewald summations. Simulations usually consisted of about 30,000 
cycles for equilibration plus a further 100,000 cycles dedicated to obtaining averages. One 
MC cycle typically consisted of NP/2 Monte Carlo moves, being the number of water 
molecules and P the number of replicas of the system. The configurational space was ex- 
plored by using four types of movement attempts: translation of one single bead of one 
molecule (30%), rotation of a single bead of one molecule (30%), translation of a whole ring 
(20%) and rotation of all the beads of a given molecule (20%). The maximum displacement 
or rotation was adjusted in each case to obtain a 40% acceptance probability. After each 
NP/2 Monte Carlo moves one volume change was also attempted. The maximum volume 
change was adjusted to obtain a 30% acceptance probability. 

V. RESULTS 

Before presenting the results for bulk water, a preliminary check was undertaken to show 
that the rotational propagator was indeed able to reproduce the rotational kinetic energy 
of a free asymmetric top, for which the rotational energy can be analytically computed via 
evaluation of the partition function. The comparison between the simulation data with the 
analytical expansion is excellent both for water and tritiated water (see Fig. [1]). It can also 
be observed that for temperatures above approximately 50K the quantum rotational energy 
of the free rotor is almost equal to the classical value (3/2A;bT), which means that for an 
isolated rigid water molecule quantum effects are only important below this temperature. 
However, in condensed matter the situation is different because there are inter-molecular 
forces, in this case hydrogen bonds, that hinder the rotation of the molecules and lead to 
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FIG. 1. Kinetic rotational energy of the isolated H2O (filled red squares) and ^H20 (open blue 
circles) molecule as a function of temperature. There is a good agreement between the path 
integral simulations and the rotational energy obtained from the theoretical partition function of 
an asymmetric top for H2O (solid line) and '^H20 (dashed line) geometry. For clarity, the rotational 
energy of ^H20 has been shifted O.lkcal/mol in the y-axis. 

the appearance of noticeable quantum effects at much higher temperatures. 

Before performing simulations for the hquid and solid phases, we need to choose the 
number of replicas, P, that are to be used in the simulations. The classical limit corresponds 
to one single replica, whereas the quantum limit is approached as the number of replicas 
tends to infinity. However, simulations can only be performed for a finite number of replicas. 
In practice the number of replicas is chosen so that it is small enough that simulations do 
not become prohibitively expensive, but high enough so as to capture the main contribution 
of the quantum effects. P represents a compromise between statistical convergence and 
theoretical accuracy, so a study of how the desired property converges with the number of 
replicas has to be carried out. We examined the convergence of the potential energy and 
the total energy as a function of the number of replicas for liquid water at T=298K and 
p=lbar (see Fig. [2]). It can be seen that both the potential and total energies increase with 
the number of replicas. There is a large increase for small number of replicas and then both 
magnitudes reach a plateau above 5 or 6 replicas. The quantum limit would be obtained by 
plotting the total energy as a function of the inverse of the number of replicas and taking 
the limit to infinite P. We found that the total energy is lower than the value at P —t- 00 by 
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FIG. 2. Convergence of the potential energy U (circles) and the total energy E (squares) as a 
function of the number of replicas P for liquid TIP4PQ/2005 water at T=298K and p =lbar. The 
dotted and dashed lines are only guides to the eye. 

about 3% for P =5 and by about 2% for P =7. In view of this we have chosen to use P =5 
replicas at room temperature. Other authors have also used a similar number of replicas 
for water at room temperaturei^i^Si^i^. For other temperatures, the number of replicas was 
chosen so as to keep the product PT approximately constant, i.e. taking P=5 at T=300K 
we arrive at PT ^ 1500K. 

A. Isotopic effects on the TMD and Cp 

One of the idiosyncratic properties of water is the existence of a maximum in density. As 
mentioned previously, the location of the TMD is affected by variations in the hydrogen mass. 
In particular, for deuterated water this maximum occurs 7K above the TMD of water and 
for tritiated water this increases to 9K. Therefore, we expect that nuclear quantum effects 
will shift the TMD of water to lower temperatures when quantum effects are implicitly 
incorporated. In addition, given that good water models reproduce the TMD of water— i^, 
it would be interesting to check whether the TIP4PQ/2005 model is also able to reproduce 
the experimental TMD. 

The equations of state of water, deuterated water, tritiated water and classical water 
were calculated at p =lbar. NpT PIMC simulations were performed at six different temper- 
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atures in each case. As the maximum in density is reflected in the third significant figure of 
the density, especially long simulations are required to reduce the statistical error, so each 
simulation consisted of at least 3 million MC cycles. The computational cost of the simula- 
tions was reduced by using the reaction field method^^, rather than Ewald summations, to 
account for the long range electrostatic forces. It has been shown that Ewald summation 
and reaction field provide very similar results for liquid water— and for Stockmayer fiuids^. 
For the system sizes studied in this work (A^=300-360) the reaction field technique yielded 
slightly higher energies and densities than the Ewald summations, but the location of the 
TMD was unchanged^. 

The results are shown in Fig|3l The location of the TMD was obtained by fitting the 
simulation data to a quadratic or cubic polynomial. The TMDs obtained are given in Table 
mil The results show that, using quantum simulations, the TIP4PQ/2005 model predicts 
that the TMD of water occurs at 284(2)K, only 7K above the experimental result. We have 
seen that upon increasing the number of replicas the TMD shifts to 280 (2) K—, indicating 
that simulation results become even closer to experimental results when more replicas are 
used. The results for deuterated water and tritiated water show that TIP4PQ/2005 is 
also able to qualitatively reproduce the shift to higher temperatures when the mass of 
the hydrogen isotope increases, in line with experimental observations. By considering the 
results for water, deuterated water and tritiated water (using P = 5 for the three molecules) 
one can see that the location of the TMD shifts 8K for deuterated water, and 12K for 
tritiated water, with respect to that of water. These are only slightly larger than the 7K 
and 9K, respectively, found experimentally. With regards to density, quantum effects affect 
differently the density depending on whether we are above or below the TMD. For high 
temperatures the number density increases with mass whereas for low temperatures the 
TMD the number density decreases when the mass is increased (see Fig. [3]). 

It is also interesting to study the shift in TMD when going from quantum water to 
classical water. Our simulations predict that the TMD can change as much as 30K when 
quantum effects are included. It seems that this is a typical shift for rigid non-polarisable 
models (a similar result was found for the TIP5P model^^). Similarly de la Pena, Razul and 
Kusalik estimated a shift in the melting point of ice Ih for the rigid TIP4P model of about 
35K when nuclear quantum effects are included^. However, the shift in the TMD between 
the quantum and classical limit could well be different for different types of potentials (for 
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example, flexible or polarisable models). In particular, it has been found that for the flexible 
polarisable TTM2.1-F model^ the explicit inclusion of quantum effects left the location of 
the TMD unchanged^. Further work is needed to clarify this. The TMD increases as the 
molecule becomes more and more classical (i.e. as the mass of the hydrogen isotope is 
increased). The number densities of water, deuterated water and tritiated water obtained 
from PI simulations along the room pressure isobar are shown in Fig. |3l Results from 
classical simulations along this isobar are also presented. As can be seen in Fig. [3] the 
number density at the maximum is hardly affected by the mass of the hydrogen isotope; 
differences between the number densities at the maximum are within the estimated error 
bar. However, it seems that the number density at the maximum first decreases slightly 
on going from water to deuterated and then tritiated water, and then increases a little in 
the classical limit. Experimentally it has been observed that the number density at the 
maximum decreases by ~ 0.7% on going from water to tritiated water (see Table [T]). PIMC 
simulations, to a lesser extent, also reflect this decrease (~ 0.2%). We have shown in previous 
work that the internal energy of water is non-linear when plotted as a function of the inverse 
of the mass of the hydrogen isotope'^^. The same also seems to be true here for the density 
of water at the maximum. This indicates that the behaviour of classical water cannot be 
obtained from a simple extrapolation of results obtained for water, deuterated water and 
tritiated water. 

Another interesting calculation is the difference between the TMD and the melting tem- 
perature obtained via classical and quantum simulations (AT = Ttmd — Tmeit)- It has been 
found that for classical simulations of rigid non-polarisable models the TMD is situated 
about 30K above the melting temperature (i.e., AT=30K)^, which is much higher than 
the AT=4K found experimentally. Quantum simulations for the flexible q-TIP4P/F model 
predict that the difference between the TMD and the melting point is also about 30K— . 
Preliminary direct coexistence simulations^!^ of the melting temperature together with the 
TMD calculations presented before indicate that the difference might be about 20-22K for 
the TIP4PQ/2005 model, which improves upon the classical prediction, but that is still 
far from the 4K found experimentally. This suggests that although the inclusion of nuclear 
quantum effects reduces the value of AT other features of real water need to be incorporated 
in the model, such as polarisability, in order to quantitatively reproduce the experimental 
difference. 
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System TIP4PQ/2005 Expt. 

H2O (P=5) 284(2) 277.13 

^HsO (P=5) 292(2) 284.33 

^HsO (P=5) 296(2) 286.55 
Classical H2O 307(2) 

TABLE III. Temperature of maximum density at p =lbar as obtained from PIMC simulations of 
the TIP4PQ/2005 model for several water isotopes (temperatures are given in kelvin). 
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FIG. 3. Isotopic effects on the TMD of water along the room pressure isobar. Number densities 
(i.e number of molecules per unit of volume) as a function of temperature at room pressure are 
presented. 



From the simulations performed along the 1 bar isobar it is straightforward to evaluate 
the heat capacity at constant pressure for water as well as other water isotopes (Cp = 
H being the enthalpy). It has been found experimentally that the heat capacity of liquid 
water is considerably affected by the isotopic substitution of the hydrogen atom (see Table 
[T]). In particular, the heat capacity is about a 10% higher for deuterated water than for 
water at room temperature^. The heat capacities obtained from PIMC simulations with 
TIP4PQ/2005 model for liquid water, deuterated water, as well as classical water simulated 
with the TIP4P/2005 model are shown in Fig. SI It can be seen that quantum simulations 
with the TIP4PQ/2005 model are able to quantitatively reproduce the heat capacity of liquid 
water for a broad range of temperatures. In addition, the increase in the heat capacity upon 
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isotopic substitution is also quite accurately captured. The agreement between simulations 
and experimental data for deuterated water is quite remarkable. Classical simulations with 
the TIP4PQ/2005 model are about 25% higher than those for quantum water using the same 
model^^, which again indicates that Cp is significantly affected by the inclusion of quantum 
effects. One might think that a good description of the heat capacity of liquid water at 
room temperature could be obtained by using a classical description with a model in which 
quantum effects are implicit though the parametrisation of the model. However, to the best 
of our knowledge, none of the rigid non-polarisable models proposed so far has been able 
to provide a quantitative description of the heat capacity of liquid water—. As an example, 
results for TIP4P/2005 are also shown in FigJH which are in poor agreement both with 
experiments and with quantum simulations. This suggests that quantum effects need to be 
incorporated explicitly in order to obtain a quantitative description of the heat capacity of 
water. This is further supported by the finding that PIMC simulations with TIP4PQ/2005 
also reproduce the experimental heat capacity of ice \h from very low temperatures up to 
room temperature (results for ice \h are shown later). These results demonstrate that the 
main contribution can be captured using a rigid model and that the contribution from 
the intra-molecular degrees of freedom is small. This is not unexpected; intra-molecular 
vibrations exhibit very high frequencies (~ 3000 cm~^) so at room temperature only the 
ground state is populated and therefore there is little or no contribution to Cp from the 
intra-molecular vibrations. 

B. Equation of state of ices 

A significant deficiency of classical simulations is the inability to reproduce the equation 
of state of solids at low temperatures. One of the consequences of the third law of thermody- 
namics is that the thermal expansion coefficient should tend to zero at zero temperature^, 
which is equivalent to saying that the density should remain constant at low temperatures. 
However, in previous work, we have seen that classical simulations of both TIP4P/2005^- 
and TIP4P/Ice^ are unable to reproduce the correct curvature of the equation of state of 
ices at low temperatures^. This was attributed to the fact that quantum effects become 
increasingly important as the temperature decreases. To see whether the description of ice 
at low temperatures could be improved when quantum contributions were incorporated, we 
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FIG. 4. Heat capacity at constant pressure for water and other water isotopes as a function of 
temperature at p=lbar calculated by means of PIMC simulations with the TIP4PQ/2005. The 
heat capacity for classical water was simulated with the TIP4P/2005. Experimental results for 
water~ (filled circles) and deuterated water®^ (filled squares) are also shown for comparison. 



performed PIMC simulations using the TIP4PQ/2005 model in order to obtain the equation 
of state of ice II. The results of quantum simulations with TIP4PQ/2005 together with the 
results of classical simulations with TIP4P/2005 and the experimental data of Fortes et al7- 
are shown in Fig. [51 These results show that a good agreement with the experimental data 
can be obtained when quantum contributions are explicitly included; the equation of state 
has the same curvature as the experimental curve, now in concordance with the third law 
of thermodynamics. The same was also found to be true for ice and for hydrate si—. 



C. Structure of ices 

It is usually found that classical simulations using simple models of water tend to over- 
estimate the height of the first peak in the oxygen-oxygen distribution function for both 
liquid water— and for ice I/i— . It is well known that quantum effects lead to less structured 
liquids^'^'^'^ and solids^ and so one might think that this can be corrected by including 
quantum effects. Radial distribution functions for liquid water and ices Ih and II obtained 
from classical simulations for TIP4P/2005 and quantum simulations with TIP4PQ/2005 at 
relatively high temperature (i.e. T >250K) and room pressure are presented in Figs. O [7] 
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FIG. 5. Equation of state of ice II at p =lbar as calculated from PIMC simulations with 
TIP4PQ/2005 and from classical MC simulations with TIP4P/2005. Experimental data are also 
shown for comparisoit^. 



and [HI Differences are visible although relatively small, being larger for the oxygen-hydrogen 
and hydrogen-hydrogen distribution functions, which is not unexpected since quantum ef- 
fects are mainly due to the hydrogen mass. The situation is different at low temperatures. 
Both classical and quantum simulations were performed for ice II at T =100K. The re- 
sults are given in Fig. |8l The first peak of the oxygen-oxygen distribution function is 
considerably lower for quantum simulations with TIP4PQ/2005 than for classical simula- 
tions with TIP4P/2005. As far as we know, as yet there are no experimental data for 
the atomic distribution function of ice II at this thermodynamic state, but it is expected 
that quantum simulations provide a better description of the structure at low temperatures 
(as was the case for ice I^— , the only ice for which the oxygen-oxygen atomic distribution 
function has been experimentally measured). The effects are much larger if one examines 
the oxygen-hydrogen and hydrogen-hydrogen distribution functions, for which differences 
between classical and quantum simulations extend further than the first peak. In partic- 
ular, the classical hydrogen-hydrogen distribution function exhibits a large number of well 
defined peaks which become considerably smoother in the quantum limit. The results for 
ice Ih and II also suggest that nuclear quantum effects might affect more significantly the 
hydrogen-hydrogen atomic distribution function in proton ordered ices, such as ice II. 
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FIG. 6. Atomic distribution function of liquid water at 298K and 1 bar as calculated from clas- 
sical MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with 
TIP4PQ/2005 (red solid line). Experimental data (black dotted line) are also showit^. 

D. Thermodynamic coefficients for ice 1^ 

In previous work, we demonstrated that classical simulations using both the TIP4P/2005 
and the TIP4P/Ice models were unable to provide a good description for many thermody- 
namic coefficients for ice Ih—- In particular, it was shown that classical simulations resulted 
in a poor description of the heat capacity at constant pressure and of the thermal expansion 
coefficient. The thermal compressibility, on the other hand, was described reasonably well. 
We checked whether the description of some of these thermodynamic coefficients could be 
improved by performing quantum PIMC simulations with the TIP4PQ/2005 model. The 
heat capacity at constant pressure and the thermal expansion coefficient can be calculated 
from the simulations that trace out the room pressure isobar (data shown in FiglS]). The 
heat capacity was obtained by fitting the enthalpy to the function H = a + bT"^ + cT^ and 
differentiating this fit with respect to the temperature. The isothermal compressibility was 
calculated by performing simulations at p =-500, -250, 0, 250 and 500 bars for temperatures 

23 




2 3 4 5 6 7 
r/A 



FIG. 7. Atomic distribution function of ice Ih at 250K and 1 bar as calculated from classi- 
cal MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with 
TIP4PQ/2005 (red sohd hue) 



between lOOK and 250K. The density along each of these isotherms could be nicely fitted to 
a straight line. The isothermal compressibility was computed by differentiating the density 



). Once the thermal expansion coef- 

T 



with respect to the pressure from the fit {kt = ^ §^ 

ficient and the isothermal compressibility are known, the pressure coefficient {Pv = ^Iv-^ 
can be readily computed via /Sy = a/i^T- 

The thermal coefficients as obtained from PIMC simulations using the TIP4PQ/2005 
model, as well as those form classical MC simulations using the TIP4P/2005 model'^^ are 
shown in Fig. O The thermal coefficients derived from the experimental equation of state 
of Feistel and Wagner are also shown for comparison^^. The results indicate that, except 
for the thermal compressibility, for which a similar accuracy is obtained in quantum and 
classical simulations, quantum simulations considerably improve the description of the ther- 
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FIG. 8. Atomic distribution function of ice II at lOOK and 250K and at 1 bar as calculated from 
classical MC simulations with TIP4P/2005 (blue dashed line) and quantum PIMC simulations with 
TIP4PQ/2005 (red soHd hne) 



modynamic coefficients. In particular, the heat capacity at constant pressure, which was 
not reproduced by classical simulations at any temperature, is now nicely reproduced from 
room temperature all the way down to zero kelvin. This is in line with the results for liquid 
water, where again the experimental heat capacity could only be reproduced by the explicit 
inclusion of quantum effects. These results strongly suggest that quantum effects are crucial 
when it comes to describing the heat capacity of either liquid water or ice Ih- 

Regarding the coefficient of thermal expansion a, quantum simulations using the TIP4PQ/2005 
model give a much better description than classical simulations with the TIP4P/2005 model. 
Even though at room temperature classical simulations predict a value closer to the experi- 
mental data, quantum simulations provide a better overall description over the whole range 
of temperatures. Importantly, the thermal expansion coefficient now tends to zero at zero 
kelvin, as it should according to the third law of thermodynamics. Finally, the description of 
the pressure coefficient {Pv) is also considerably improved when including quantum effects. 
This simply reflects that the thermal expansion coefficient is improved in a quantum de- 
scription of the system, since /3v = a/ ^t- The isothermal compressibility is little affected by 
quantum effects; almost all the change in the pressure coefficient is due to a good description 
of the thermal expansion coefficient. 
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FIG. 9. Thermodynamic coefficients {kt, a, /3y and Cp) of ice at p=lbar, as calculated from 
classical simulations with TIP4P/2005 and from quantum simulations with TIP4PQ/2005. Exper- 
imental data are also shown for comparisort^. 



E. Relative energies of ices at zero kelvin 

Finally we have also computed the relative energies between various ice phases at zero 
temperature. It has been found that some classical water models result in a rather good de- 
scription of the phase diagram of water— i^"^. However, there is still room for improvement. 
For example, it has been found that usually ice II is over-stabilised with respect to ice Ih, 
for some models, so much so that ice II completely removes ice Ih from the phase diagrairt^. 
A preliminary outline of the phase diagram for a particular model can be obtained by esti- 
mating the coexistence pressures between the competing solid phases at zero temperature. 
At zero kelvin phase transitions occur with zero enthalpy change, so a calculation of entropy 
is avoided. Assuming that the change in energy and density between two solid phases is 
almost independent of pressure at zero kelvin (which is indeed a rather good approximation 
for ices), the calculation of coexistence pressures between two ices at zero temperature can 
be estimated from^-: 
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(50) 

p=0 



Therefore, by simply calculating the energy and density of the solid phases at zero temper- 
ature and zero pressure one can obtain a reasonable estimate of the coexistence pressure at 
zero temperature. 

The properties at zero temperature were computed for ices Ih, II, III, V and VI. Empty 
hydrates structures si, sll and sH, which have been shown to be the stable solid phases at 
negative pressures^, were also considered. Simulations were performed along the zero bar 
isobar in the temperature range from 250K to 77-lOOK. The energy at zero temperature 
was obtained by fitting the data to the function E = a + hT"^ + cT^, from which one can 
estimate E{T = OK). The energies obtained using this procedure are represented in Fig. 
[TOl Energies are given relative to the energy of ice Ih, which experimentally is the most 
stable phase at zero temperature and at zero pressure. The results show that both classical 
MC simulations using TIP4P/2005 and quantum PIMC simulations using TIP4PQ/2005 
predict that ice Ih is the most stable phase, in agreement with experimental results. It is 
also observed that the relative energies of ices II, III, V and VI change when quantum effects 
are explicitly taken into account. In particular, ice Ih is destabilised with respect to ice II by 
about 0.2kcal/mol, so that now the relative stability of ice II with respect to ice Ih is much 
closer to the experimental value. This indicates that quantum effects are larger in ice Ih than 
in ice II, resulting in a de- stabilisation of the former. Ices III, V, and VI are also stabilised 
with respect to ice Ih, but to a smaller extent (by about 0.1 kcal/mol). Finally, the relative 
stability of the empty hydrate structures are not appreciably changed. Taking everything 
into account, we can identify three different families of ices according to the importance of 
quantum effects. The first family includes ice Ih and the empty hydrate structures si, sll 
and sll, which are influenced the most by quantum effects. The second family will be that 
formed for ices III, V, and VI, and finally, ice II, which is the least affected by quantum 
effects, forms the third family. 

The reason why quantum effects make distinct contributions to the various ice phases can 
be understood by looking at the geometrical arrangement of the four molecules that form a 
hydrogen bond with a central one. These molecules form a nearly perfect tetrahedron in ice 
Ih and a slightly deformed tetrahedron with deviations of about 10 degrees in the hydrates. 
On the other hand the deviations from the perfect tetrahedron are of about 30 degrees for 
the ices II, III, V and VI. As a result, the strength of the hydrogen bond is larger in ice Ih 
and the hydrates than in the remaining ices and, as can be seen in Eq. HTJ this results in 
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FIG. 10. Estimate of the relative energies at zero temperature of the ice phases as calculated from 
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an increased value for the average forces and torques, boosting the quantum influences. 

Once the energies at zero temperature have been calculated, the coexistence pressures 
can also be obtained. The results are given in Table IIVI As can be seen, the coexistence 
pressure between ices Ih and II is the most affected coexistence line. It decreases from about 
2090 bar in classical simulations (TIP4P/2005) to about 195 bar in quantum simulations 
(TIP4PQ/2005), resulting in a much better agreement with the experimental results. The 
coexistence lines between the remaining ices (II, III, V and VI) are also affected, but to 
a lesser extent. In general, quantum results are closer to the experimental data than the 
results obtained using classical simulations. Finally, coexistence pressures between Ih and 
the hydrates or between hydrates themselves are largely unaffected by quantum effects. 

In summary, phase transitions between solid phases belonging to the different families are 
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Phase TIP4P/2005 TIP4PQ/2005 Expt. 





2090 


195 


140(200) 


ih-ni 


3630 


2727 


2400(100) 


II-V 


11230 


15731 


18500(4000) 


II-VI 


8530 


10935 


10500(1000) 


III-V 


3060 


1998 


3000(100) 


V-VI 


6210 


6848 


6200(200) 


h-sl 


-4174 


-3948 




Ih-sII 


-3379 


-3249 




Ih-sH 


-4072 


-3933 




sll-sl 


2787 


2267 




sII-sH 


-7775 


-7557 





TABLE IV. Coexistence pressures (in bar) at zero temperature obtained from quantum PIMC 
simulations with TIP4PQ/2005 and classical MC simulations with TIP4P/2005. These data are 



taken from Refs. 
comparison. 
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71 



76l . and 



78l . Experimental data taken from Ref. 177) is also given for 



the most affected by the inclusion of quantum effects, whereas phase transitions involving 
phases of the same family are in general less affected. Therefore, the explicit inclusion of 
quantum effects is crucial if one wishes to reproduce the phase transitions between solid 
phases belonging to different ice families, especially the transition Ih-ll, whereas usually 
phase transitions between ice phases of the same family can be calculated by means of 
classical simulations in conjunction with a good classical model. 



VI. CONCLUSIONS 

In this work it has been shown in detail how the formulation of the path integral for rigid 
rotors, derived by Miiser and Berne^^, can be applied to water. Using this formulation, a 
large number of properties of liquid water, ices and hydrates have been studied by PIMC 
simulations using the TIP4PQ/2005 model, which was specifically designed to use within 
quantum simulations^. For liquid water, isotopic effects on the TMD and on the heat 
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capacity have also been considered. The results show that in general a better description 
of water is obtained when quantum effects are included, although some properties can also 
be reasonably described in classical simulations with a good classical model. In addition, 
PIMC simulations with TIP4PQ/2005 reproduce the experimental isotopic effects on the 
TMD and the heat capacity of liquid water. 

Quantum effects have been found to be crucial when it comes to reproducing many prop- 
erties of water and ices. In particular, quantum effects have been found to be most important 
with regards to the properties of ices at low temperatures, which is not entirely unexpected. 
In particular, classical simulations fail to reproduce the curvature of the equation of state at 
low temperatures found experimentally and imposed by the third law of thermodynamics^. 
Our results show that this can be corrected by including quantum effects resulting in physi- 
cally agreeable equations of state for ices Ih— and II over a quite broad range of temperatures. 
Classical simulations overestimate the first peak in the oxygen-oxygen radial distribution 
function of ice Ih at 77K, which again is brought into to agreement with experiment when 
quantum simulations are performed^. The same behaviour was found for ice II and is ex- 
pected to occur for other ices. As a result of the better description of the equation of state 
the coefficient of thermal expansion at low temperatures of ice Ih is also improved when 
quantum effects are included. In addition to the properties of ices at low temperatures, 
quantum effects have also been found to be important when it comes to reproducing the 
heat capacity of ice Ih and water at all temperatures. In addition, isotopic effects on the 
heat capacity of liquid water have also been captured. 

We also found that the magnitude of quantum effects is different for different ices and, 
therefore, they need to be included if one wishes to improve the description of phase transi- 
tions. In particular, it has been found that ices can be classified into three different families, 
according to the importance of quantum effects: the first family is formed by ice Ih and 
the hydrates structures si, sll and sH, for which quantum effects are the largest, the second 
family comprises ices III, V and VI, and the third family is formed by ice II, for which quan- 
tum effects are the smallest. Phase transitions between ices belonging to different families 
change when quantum effects are included, whereas transitions between ices belonging to 
the same family are only slightly affected by quantum effects. As quantum effects are also 
different for liquid water and ice Ih, the melting point of water is also affected by quantum 
effects; it has been found in previous works that the melting point shifts to lower tempera- 
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tures when quantum effects are included^S'^. In general, quantum affects sfiould also affect 
any property that involves two phases for which quantum effects are different. For example, 
in previous work it has been shown that quantum effects improve the description of the 
enthalpy of vaporisation^^ and the sublimation enthalpy^. 

However, it has been found that other room temperature properties, although also af- 
fected by quantum effects, can be properly described using a classical model. This can be 
explained because classical models are usually fitted to reproduce some experimental data 
at room temperature, so in some way quantum effects at this temperature are implicit in 
the model. For example, it has been found that the structure of liquid water and ices Ih and 
II above 250K is reproduced with similar accuracy in classical and quantum simulations (al- 
though only quantum simulations can reproduce isotopic effects)^. The densities of ices can 
also be reproduced at room temperature with good accuracy using classical simulations^. 
The isothermal compressibility of ice Ih shows little improvement with respect to classical 
simulations at all temperatures, which also might be related to the fact that quantum effects 
are influenced little by pressure. 

With regards to the TMD, several classical models were proposed that reproduce the 
experimental TMD^*^, although the temperature difference between the TMD and the 
melting temperature is largely overestimated (they usually predict a 30K differenced^ in- 
stead of the 4K found experimentally). Preliminary PIMC calculations indicate that this 
difference might be reduced for TIP4PQ/2005 to ~ 22K, which means that, although some 
improvement is achieved, other features of real water need to be included, such as polaris- 
ability and flexibility, in order to obtain a better agreement with experiment. 

By using a rigid model, we have ignored the influence of quantum effects in the intra- 
molecular degrees of freedom. Despite this seemingly drastic approximation the results 
presented here seem to indicate that for many properties the main quantum contributions 
arise from the inter- molecular degrees of freedom. Comparison between quantum simulations 
of rigid and flexible models will be be very useful to quantify the relative importance of 
quantum effects on inter-molecular and intra-molecular degrees of freedom. 
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FIG. 11. Reference system fixed in the water molecule adopted in this work. This convention is 
usually denoted as the bca convention, because the x,y and z axis coincide with the b, c and a 
principal axes of inertia (thus la < h < Ic)- The origin of the coordinate system is located at the 
centre of mass of the water molecule. 



Appendix A. Asymmetric top eigenfunctions 

The rotation of an asymmetric free rotor can be obtained by solving the Schrodinger 
equation: 

where T''°* is the Hamihonian associated with the angular momentum: 

^™* = ^ + ^ + ^ (52) 

In this equation L^., Ly and are the three components of the angular momentum and l^x-, 
lyy and Izz are the three components of the momentum of inertia. To solve this equation 
it is convenient to choose a reference system so that the x, y and z axis are located along 
the three principal axes of inertia, denoted as a, h and c. We adopt the convention that 



la < Ih < Ic- Note that there is no unique way to identify x, z with a, 6, c (see Refs. 



20 



andlSd). For example, one could associate x with a, y with 6, and z with c, which is usually 
referred to as abc convention. Alternatively we could choose to identify x with 6, y with c 
and z with a, which is usually known as hca convention, and it is this convention that was 
used in this work (see Fig. [TTj) . The choice of axis is highly relevant since it defines the 
Euler angles that appear in the three components of the angular momentum. 
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In an asymmetric top all three moments of inertia are distinct {la ^ h ^ Ic)- In this 
situation the Hamiltonian commutes with and with L^, but not with Lc'. 

[P, H'"'] = [L„ R'"'] = [Le, H'"'] ^ (53) 

Therefore, the eigenfunctions of the Hamiltonian will also be eigenfunctions of and L^, 
but not of Lc'. 

= J(J+ 1);^^^ J = 0, l,...,cx) (54) 

L^^ = M;i^ M = -J,...,0,..., J (55) 

The solutions for the Schrodinger equation for an asymmetric top will be denoted as \JMK). 
The integer K is not a true quantum number (i.e., it does not quantise any observable) it is 
simply a number used to label the (2 J + 1) possible values of the energy available for each 
value of J and M. The functions | JMK) can be obtained by expanding them in a basis set 
formed by the eigenfunctions of the spherical top {\JMK)): 

\JMk)=Y,A'i'^\JMK) (56) 

K 

Using the hca convention, the energies E^^^^'' of the asymmetric top and the coefficients ^^^^ 
(i.e, the eigenvectors A^*^^) can be obtained solving the following secular determinant (see 



Refs. 



79 



and 



80l ) for each value of J and M : 
1 



Hkk = ^{B + C) [J{J + 1) - K^] + AK^ 
Hkk±2 = - C)[J{J + 1) - K{K ± l)fV{J +1)-{K±1){K± 2)f' (57) 

where K ranges from —J to +J. The remaining elements of the determinant are zero. A, 
B and C are the rotational constants, A = -Ar, B = and C = -ri-. Note that since 
la ^ h ^ I CI it follows that A > B > C . This determinant has dimensions of (2 x J + 1). 
Therefore, (2 x J + 1) eigenvalues are obtained for each value of J and M, which are, in 
general, all different. These (2 x J + 1) energy levels are labelled with the K index. However, 
as M does not appear in the determinant, there is a (2 x J + 1) degeneracy in the energy 
associated with M. Note that if we had chosen a different convention for the reference 
system, e.g. the abc, the secular equation that we would need to solve would be different 



(see, e.g., Refs. 



19|, 



79 



, and 
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Appendix B. Wigner functions 



The Wigner functions are given hy^: 

diM = [{J + M)\{J - M)!(J + K)\{J - K)\f' (58) 

^ iJ-M~x)KJ + K- x)Kx)Kx + M-K)\ 

where the sum over x is restricted to those values that do not lead to negative factorials. 

Appendix C. Obtaining the Euler angles of a water molecule from the Carte- 
sian coordinates of its interaction sites 



In our implementation of the MC algorithm we used Cartesian coordinates. This is 
convenient when it comes to computing the potential energy between two water molecules. 
However, to evaluate the density function of the quantum rotational energy, the orientation 
of the molecules must be specified in terms of Euler angles. Therefore, we need a procedure 
to obtain the Euler angles that define the orientation of a given molecule from the Cartesian 



81 



e 



coordinates of its interaction sites. We used Euler angles (6', 0, x) as defined in Ref 
varies from to tt and (j) and x go from to 27r. Let us denote X, Y, Z as the orthogonal axes 
of a laboratory frame fixed in the space. Let us assume that we associated three orthogonal 
axes to the molecule, namely x, z, which define the body frame (with its origin located 
at the centre of mass of the molecule). We shall assume that both sets of orthogonal axes 
are right handed. The orientation of the molecule can be defined by three Euler angles. 
The three Euler angles are defined by the operations required to move the molecule from 
an initial orientation, where x, y, z are coincident with X, F, X (having both set of axes a 
common origin) to its current configuration (where the two set of axes also have a common 
origin). Rotations are counterclockwise. First, a rotation is performed about the Z-axis by 
an angle 0, so that the axes X, F, Z change to x\ y\ Z. Secondly, we perform a rotation 
by an angle 9 about the y'-axis obtained from the previous rotation {x',y',Z — )• x",y',z'). 
Finally, a new rotation by an angle x is performed around the 2;'-axis of the frame fixed in 
the body ( x", y', z' to x, y, z). Using this convention, the coordinates in laboratory frame of 
a site of the molecule, X, Y, Z (i.e., R ) can be obtained easily from the coordinates of that 
site in body frame, Xh,jh,Zh (i.e., rb ) using the expression : 
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X 
Y 



COS COS d COS X ~ sin sin x — cos cos sin x — sin cos x cos sin Q 
sin cos 6* cos x + cos sin x ~ sin cos Q sin x + cos cos x sin sin 
^ — sin 6' cos X sin 6' sin x cos y 



(59) 



Let us assign a particular body frame to the molecule of water. Taking the hca convention, 
the body frame of the water molecule is chosen so that the h principle axis (i.e., that 
associated with Jfc) is assigned to x, c is assigned y and a is assigned to z (remember that 
la h Ic)- Under this convention, the b axis lies along the H-O-H bisector, the c axis is 
perpendicular to the molecular plane and the a axis is parallel to the line connecting the two 
hydrogen atoms. With this choice the coordinates of the oxygen are (a, 0, 0) (the centre of 
mass in water is located upon the H-O-H bisector, slightly below the oxygen atom, and we 
shall define the x direction such that a is positive). The coordinates of the hydrogen atoms 
are be given by (/3, 0, —7) and (/3, 0, 7). Notice that the value of (3 should be negative, since 
a was taken to be positive. We denote {Xo,Yo, Zq), {Xh^^Yh^, Zh^) and (X//2, Y^g) -^^2) 
as the coordinates of the oxygen and of the hydrogens, respectively, in the laboratory frame, 
whose origin is fixed to be the centre of mass of the molecule. As mentioned before, the 
centre of mass of the molecule coincides with the origin of the body frame fixed in the 
molecule. For the oxygen, the coordinates in the body frame fixed in the molecule (a, 0, 0) 
are related to those of the laboratory frame through Eq. [591 The multiplication of matrices 
in Eq. |59] leads to three equations: 

Xo = a(cos0cos6'cosx — sin0sinx) (60) 
Yo = a(sin0 cos 6* cos X + cos sin x) (61) 
Zo = —a sin 6' cos x (62) 
Analogously, the coordinates of the hydrogens in the laboratory frame are given by: 

= P (cos cos 6 cos X — sin sin x) ~ 1 cos sin 6 (63) 

= P (sin cos 6 cos x + cos sin x) ~ 7 sin sin 6 (64) 

Zh^ = —f3 sin 9 cos x ~ 7 cos 9 (65) 

Xh2 = P (cos cos 9 cos X — sin sin x) + 7 cos sin 9 (66) 
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= (sin cos 9 cos x + cos sin x) + 7 sin sin ^ (67) 

= — /3 sin 6' cos x + 7 cos 6* (68) 

In summary, we have obtained nine equations to determine the three Euler angles (0 and 
X go from to 27r and, therefore, their value can only be unambiguously obtained from 
the knowledge of both their sine and cosine; whereas 6 varies from to vr so that it is only 
necessary to know its cosine). The Euler angles for an instantaneous configuration where 
the atomic coordinates are {Xq,Yq, Zq), {Xh^.Yh^, Z^^) and (X^^-j , ; ) can therefore 
be obtained by solving the set of equations defined by Eqs. [00] - EHl Subtracting Eq. 
from [65] we obtain: 

cos 6 

Subtracting EqsJ66] from [63] 

cos = 

and subtracting Eqs. [03 from 



Zh2 — Z 



HI 



27 



27 sin 9 
27 sin 9 

Finally, the Euler angle x can be obtained from Eqj62] 

-Zo 



(69) 
(70) 



sm I 



cosx 



a sin 9 



and adding Eql03]and Eqj66] 



sin X = ( cos cos 9 cos x 



or, alternatively, adding Eqj6l]and Eq. [67] 

sin X = I — sin (j) cos 9 cos x + 



2/3 



sm I 



2/3 



cos ( 



(71) 



(72) 



(73) 



(74) 



In the special case that ^ = or = tt the expressions given above are not valid because the 
denominator vanishes, resulting in a singularity. The probability of obtaining exactly 6' = 
or = vr is very small during a simulation. In these special cases, alternative expressions 
can be obtained by evaluating of the rotation matrix for the particular value of 9. When 
^ = 0, sin^ = and cos 6* = 1 and, therefore, the rotation matrix (Eq. [59]) becomes 

'^cos(0 + x) — sin(0 + x) 0^ 
sin(</) + x) cos(0 + x) (75) 
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i.e, in this case, the rotation can be seen as a simple rotation about the 2;-axis by an angle 
(j)' = (j) + X- III the case 6' = there is no a unique way of defining and x individually, 
as any combination of and x having the same value of 0' (0' = + x) will provide the 
same final configuration. Here for the particular case ^ = we decided to assign = 0, and 
with this choice the sine and cosine of the x angle can be readily obtained using a procedure 
similar to that outlined above for a general case but using the rotation matrix given in Eq. 
[751 Using this procedure we obtain : 



and 



sin(x) 



cos(x) 



X 



Hi 



/3 



(76) 



(77) 



Finally, when 6 = tt, sm6 = and cos^^ = —1 so that the rotation matrix is now: 

^- cos(x - 0) sin(x - 0) ^ 
sin(x - 0) cos(x - 0) 



v 











(78) 



-V 



In this case, this rotation can be seen as a simple rotation about the 2;-axis by an angle 
0' = X — 0. Again it is not possible to assign in a unique way values of x and 0. For this 
reason we arbitrarily assigned in this case = so that x is obtained as: 

yni 



smx 



cosx 



'1^ 



(79) 



(80) 



Appendix D. Obtaining the relative Euler angles of replica t + 1 with respect 
to those of replica t. 



Let us focus on two replicas t and t + 1 of a certain molecule. Suppose that R* are the 
laboratory frame instantaneous coordinates of a certain site of replica t with respect to its 
centre of mass (the centre of mass of replica t). These coordinates can be obtained from the 
body frame coordinates of that site (r;,) using the rotation matrix Mf. 

R* = M^r, (81) 
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Inverting the previous equation one obtains 



r;. = Mr'R* (82) 

The relative Euler angles of replica t + 1 with respect to those of replica t can be computed 
by expressing the instantaneous coordinates of replica t + 1 (R*"*"^) in the reference system 
of replica t by: 

R'*+^ = Mr'R*+' (83) 

In other words, replica t+1 is rotated using the rotation matrix M^"^ to obtain its orientation 
with respect to that of replica t. The relative orientation of molecule t + 1 with respect to 
that of molecule t is given by the atomic coordinates R'*"^^. The Euler angles associated 
with this orientation, i.e., the Euler angles of replica t + 1 with respect to those of replica t, 
can be computed using the procedure described in Appendix C. 
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